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DESCRIPTION  OF  PROBLEMS 


In  a  wide  variety  of  Air  Force  applications,  highly  detailed 
information  about  atoms,  molecules,  and  their  interactions  is  required^.1-3) 

This  information  is  necessary  in  problems  ranging  from  chemical  laser 
development,  to  the  detection  and  identification  of  rocket  plumes,  to 

(1-31 

metal  clustering  and  aerosol  formations,  and  even  to  nuclear  weapons  effects. 

Probably  the  most  crucial  component  needed  to  understand  molecular 
reactions  is  the  potential  energy  surfaces  that  serve  to  describe  the 
attractions  among  the  atoms  and  mol eculesC1  ^ However,  such  information  is 
not  easy  to  obtain.  A  certain  amount  of  information  about  the  molecular 
forces  near  equilibrium  in  a  bound  molecule  is  available  from  spectroscopy. 

Some  information  about  the  potential  energy  surface  even  in  the  absence  of 
binding  can  be  provided  from  crossed  molecular-beam  experiments.  But,  in 
general,  potential  energy  surfaces  are  not  amendable  to  experimental 
determination.  Instead,  other  types  of  experimental  observations  such  as 
kinetics  experiments,  coupled  with  very  simple  theoretical  models  for  a 
surface,  are  used  to  infer  pieces  of  information  about  the  parameters  of  the 
model  such  as  what  the  activation  barrier  might  be. 

The  most  direct  approach  to  obtaining  detailed  information  about 
a  potential  energy  surface  is  offered  by  predictive,  ab  initio  quantum 
mechanical  calculations.  However,  to  make  it  feasible  to  calculate  accurate 
energy  surfaces  for  molecules,  much  better  and  more  computationally  efficient 
methods  must  still  be  developed. 

(4-151 

One  such  approach,  namely  many-body  perturbation  theory  (MBPT)  ' 
and  its  infinite-order  extensions  termed  coupled-cluster  methods  (CCM)Ol ,16-20) 
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offer  a  number  of  attractive  features  that  the  more  traditional  configuration 

(21 1 

interaction  approaches  lack.  '  During  the  first  two  years  of  this  grant 

very  efficient  computer  codes  to  perform  MBPT/CCM  calculations  were  written 

and  employed  for  the  first  time  in  large-scale  ab  initio  calculations  of 

potential  energy  surfaces  The  successes  of  this  effort  have  been  substantial. 

These  include  the  determination  of  a  complete  force-field  for  the  H20 

molecule,  including  all  force-constants  through  fourth-order,  that  is 

sufficiently  accurate  that  once  improved  experiments  were  carried  out  after  our 

calculations,  many  of  the  previously  accepted  values  fcr  the  force  constants 

(221 

were  revised  to  be  more  consistent  with  our  predictions.'  '  Also,  a  study  of 
the  binding  energies  of  the  molecules  B2Hg->-2BH3,  H3BNH3->-BH3+NH3,  and 
H3BC0l+BH3+C0  was  made  that  predict  these  binding  energies  to  within 
1  kcal/mole  of  the  accepted  experiments  for  diborane  and  borane  carbonyl, 

(14l 

and  made  a  prediction  in  the  case  of  borazane  in  the  absence  of  an  experiment.'  ' 

Earlier  experiments  which  gave  much  higher  values  for  the  binding  energies 

of  diborane  and  borane  carbonyl  than  we  computed  are  now  completely  discounted. 

Similar  successes  with  studies  of  the  isomerization  energy  and  activation 

barrier  of  HNC-*-HCNp^and  CH^NC+CH^Np^the  photodissociation  of  forma IdehydeP^ 

( 26 1 

and  various  studies  of  methanol,  methoxy,  and  the  formyl  radical'  'attest  to 

the  reliability  of  our  MBPT/CCM  methods. 

Building  upon  this  work  supported  by  the  AFOSR  we  have  now 

carried  out  extensive  studies  of  the  potential  energy  surface  for  the 

3  3 

two  inelastic  collisions,  0(  P)+H20  and  0(  P )+C02 ,  under  contract  to  the 

Air  Force  Rocket  Propulsion  Laboratory,  for  the  purpose  of  obtaining 

vibrational  excitation  cross-sections  that  are  needed  in  actual  detection 
(271 

devices.'  ' 


Despite  the  many  successes  we  have  had,  there  are  still  categories 
of  problems  that  cannot  yet  be  attacked  by  MBPT/CCM.  These  include  studies 
of  most  excited  states,  reactions  that  break  multiple  bonds,  and  applications 
to  various  kinds  of  open-shell  molecules^To  satisfy  these  additional 
requirements  it  is  necessary  to  simultaneously  develop  the  formal  theory, 
write  additional  computer  programs,  and  continue  to  make  landmark  applications 
of  our  developing  quantum  mechanical  technology.  Although  in  many  cases  the 
formal  theory  is  less  dramatic  than  the  applications,  the  continual  extension 
of  the  theory  has  a  greater  impact  on  our  ability  to  calculate  accurate 
energy  surfaces  for  whatever  categories  of  problems  might  emerge  from  the 
needs  of  the  Air  Force. 

Consistent  with  this  objective,  much  of  our  work  this  past  year 
has  been  devoted  to  formal  theory.  This  includes  optimization  of  orbitals 
within  the  coupled-cluster  framework  and  developing  additional  mathematical 
techniques  to  efficiently  solve  the  non-linear  coupled  cluster  equations. 
Additional  applications  to  a  variety  of  problems  have  also  been  accomplished. 

In  the  following,  Section  II  discusses  the  research  objectives 
of  this  grant,  and  summarizes  some  of  the  notable  accomplishments  made  in 
the  past  year.  For  the  previous  year's  effort,  we  refer  to  last  year's 
annual  report.  Section  III  lists  the  publications  and  presentations  which 
have  been  supported  by  this  grant.  Section  IV  discusses  in  some  detail 
the  idea  of  orbital  optimization  in  coupled-cluster  theory.  This  is  a 
radically  new  approach  of  substantial  scientific  interest.  Appendices  A 
and  8  are  two  manuscripts  recently  accepted  for  publication.  The  first 
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REVIEW  OF  RESEARCH  ACCOMPLISHMENTS 


The  overall  objectives  of  this  research  program  include  the 
following: 

(1)  Develop  new,  rore  accurate  and  more  efficient  ab  initio 
quantum  mechanical  methods  based  upon  MBPT  and  CCM  for 
determining  molecular  properties  and  particularly, 
potential  energy  surfaces  for  molecular  interactions. 

(2)  Implement  these  methods  in  highly  efficient,  transportable 
computer  codes,  to  enable  computations  on  potential 
energy  surfaces  to  be  made  on  an  almost  routine  basis. 

(3)  Apply  these  techniques  to  a  variety  of  problems  that 
are  of  interest  to  AFOSR,  and  that  serve  to  establish 
the  range  of  accuracy  for  MBPT  and  CCM  methods. 

In  line  with  these  overall  objectives,  a  number  of  accomplishments 
have  been  made  so  far  in  this  program.  The  accomplishments  from  the 
previous  year  are  listed  in  the  Annual  Report  for  1979.  Hence,  we  will 
summarize  only  the  additional  achievements  that  have  been  made  in  the 
past  year. 

The  main  focus  for  our  effort  this  last  year  has  been  formal, 
being  directed  toward  the  generalization  of  the  MBPT  and  CCM  theory.  In 
particular,  the  inclusion  of  monoexcited  clusters  and  the  theory  required 
for  orbital  optimization  have  been  developed  within  the  CCM  model.  Additional 
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formal  work  has  been  directed  at  the  multi  reference  MBPT/CCM  theory.  A 
series  of  applications  to  a  wide  variety  of  problems  using  our  previously 
established  computer  codes  have  also  been  made  to  demonstrate  the  applicability 
of  our  developing  methods.  A  summary  of  achievements  follows. 

A.  The  coupled-cluster  theory  and  programs  have  been 
generalized  to  include  monoexcited  clusters  (i.e.,  T-j). 

This  is  found  to  be  important  in  obtaining  correct 
potential  energy  curves  for  cases  where  a  single¬ 
determinant  reference  function  is  not  entirely 
appropriate. 

B.  For  the  first  time,  the  effect  of  optimizing  the 
molecular  orbitals  in  the  coupled-cluster  theory  has 
been  studied.  This  model  is  similar  to  a  multi¬ 
configuration  self-consistent  field  (MCSCF)  approach, 
except,  via  CCM,  one  employs  all  single,  double,  and 
quadruple  type  Cl  excitations.  Nothing  of  this  magnitude 
has  been  attempted  previously.  We  are  using  this 
method  to  study  some  unusually  difficult  molecular 
problems . 

C.  A  new  numerical  techniquei^similar  to  the  reduced 

partitioning  procedure  developed  for  eigenvalue 

( 29) 

problems  by  the  author  some  years  ago,  has  been 
generalized  to  apply  to  the  nonlinear  coupled-cluster 
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equations.  This  has  the  effect  of  greatly  increasing 
the  rate  of  convergence  when  solving  the  equations. 

This  has  enabled  us  to  obtain  accurate  results  for 
some  pathological  cases  that  could  not  have  been  obtained 
by  the  standard  iterative  approach.  (See  Appendix  B). 

D.  The  decomposition  of  formaldehyde,  ^CO,  to  radical  and 
molecular  products  and  its  rearrangement  to  hydroxycarbene 
has  ~een  studied i^Vhis  problem  is  of  substantial  experi¬ 
mental  interest  because  of  formaldehyde's  prominence  in 
combustion/plume  processes.  The  activation  barriers  and 
heats  of  reactions  have  been  obtained.  In  the  latter  case, 
agreement  with  experiment  is  within  +  2  kcal/mole.  The 
activation  barrier  predictions  support  the  Cl  results  of 
Goddard  and  Schaefer^ 30^  that  would  suggest  a  tunneling 
mechanism  for  h^CCW^+CO. 

E.  A  series  of  detailed  comparisons  of  various  MBPT  models  with 

CCD  for  the  C,  N,  and  0  atoms  and  the  H^O,  NH3,  and  CH4 

d21 ) 

molecules  were  made  this  yearr  'These  comparisons,  plus  a 
number  of  others  we  have  made,  suggest  that  the  infinite-order 
CCD  results  differ  insignificantly  from  the  fourth-order  model, 
DQ-MBPT(4),  for  most  normal  cases^^This  supports 
the  predictions  of  the  less  expensive  fourth-order  model 
for  larger  molecules. 
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F.  An  experimentalist  measured  the  isomerization  energy  of 

HNOHCN  to  be  10+1  kcal/mole.^^  This  disagreed  with  our 

(241 

theoretical  prediction  of  15+2  kcal/mole.'  '  To  attempt 
to  resolve  this  discrepancy,  we  performed  calculations  on 
this  system  for  a  series  of  approximations  and  basis  sets, 

concluding,  indeed,  that  the  value  15+2  kcal/mole  is 

f  23 1 

accurate.  '  We  believe  this  will  prove  to  be  another 
case  where  theory  has  demonstrated  that  the  experimental  value 
is  in  error. 

G.  In  addition  to  the  HNOHCN  rearrangement,  the  interesting 

(23) 

systems  LiNO-LiCN  and  BNOBCN  were  also  studied.  Essentially 
no  barrier  to  rearrangement  is  found  for  LiNOLiCN,  while 
BNC  is  found  to  be  more  stable  than  BCN.  Our  calculations 
made  predictions  of  the  thermochemistry  parameters  for 
these  molecules  which  we  hope  will  stimulate  some 
experimental  work  . 

H.  The  first  all-electron  ab  initio  coupled-cluster  and  MBPT 

calculations  of  benzene  were  made  this  year.  This  work 

demonstrates  that  a  molecule  of  this  size  has  at  least  a 

20  percent  error  in  its  correlation  energy  due  to  the 

( 32) 

neglect  of  Cl  type  quadruple  excitations.'  '  This 
emphasizes  the  importance  of  using  methods  like  MBPT/CCM 
that  properly  include  such  higher  order  excitation  effects 
if  reliable  quantum  mechanical  calculations  are  to  be 
possible  for  larger  molecules.  (See  Appendix  A.) 


-1 
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I.  Additional  calculations  at  the  MBPT/CCM  level  on  a 
variety  of  systems  including  Li^,  CH^O,  HNO,  and 
other  unusual  molecules  are  also  being  made  to 
predict  structures,  thermochemistry,  and  other 
properties  like  Jahn-Teller  distortions. 
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IV. _ SYNOPSIS  OF  ORBITAL  OPTIMIZATION  IN  COUPLED  CLUSTER  THEORY 

A  new  idea  which  has  been  developed  over  the  second  year  of  this 
grant  is  the  optimization  of  molecular  orbitals  while  simultaneously  carrying 
out  coupled  cluster  calculations.  This  approach  is  similar  in  philosophy 
to  multi-configuration  self-consistent  field  (MCSCF)  theory,  but  differs 
substantially  in  numerous  other  respects. 

The  traditional  MCSCF  approach  takes  a  small  number  (typically 
5-20)  configurations  which  we  will  designate  as  {D^},  composed  of  a  set 
of  molecular  orbitals ,{x. }.  The  wavefunction  then  takes  the  form 

%SCF  =  £  CkDk 

where  each  is  a  determinant  (or  symmetry  adapted  combination  of 
determinants)  of  the  general  form 

Dk  »  A(Xl(l)x2(2)...xa(i)...xb(j)...xn(n)) 

Various  replacements  (like  Xg(i)  and  Xb(j))  of  molecular  orbitals 
which  are  occupied  in  D-|  (often  the  SCF  determinant)  give  rise  to  the 
usual  Cl  single,  double,  etc.  excitations.  The  linear  coefficients  Ck 
are  optimized  by  using  the  variational  principle  in  a  small  Cl  eigenvalue 
problem.  This  is  the  multi -configuration  step.  The  SCF  part  occurs  when 
the  coefficients  { cu >  in  the  molecular  orbitals  (x^ > » 

xi  =  JVu 

for  $u  some  primitive,  usually  atomic  orbital  basis  set,  are  also  simultaneously 
optimized.  This  procedure  is  particularly  useful  in  types  of  open-shell 
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problems,  and  in  descriptions  of  bond  breaking  when  the  form  of  the 
molecular  orbitals  can  change  significantly. 

However,  there  are  a  number  of  disadvantages  with  the  MCSCF 
method.  The  small  number  of  configurations  have  to  be  carefully  chosen 
if  correct  answers  are  to  be  obtained  and  this  requires  a  difficult  and 
time  consuming  trial  and  error  procedure.  As  a  consequence,  MCSCF  methods 
have  notoriously  bad  convergence  properties. 

The  most  important  weakness  from  our  viewpoint  pertains  to  MCSCF's 
failure  to  be  size-extensive.  That  is,  the  calculations  do  not  scale 
properly  with  molecular  size.  This  occurs  due  to  the  method  using  a 
variational  Cl  step,  which  makes  it  suffer  from  the  same  weakness  as  any 
truncated  Cl,  inspite  of  optimizing  the  orbitals.  MCSCF's  failure  to  be 
size-extensive  causes  problems  ranging  from  decreasing  its  applicability 
to  large  molecules  to  correct  predictions  of  dissociation  energies.  The 
various  difficulties  encountered  with  nonsize-extensive  methods  are  documented 
in  detail  elsewhere.  ^ 

Recently,  a  series  of  papers  have  appeared  that  develop  an  approach 
for  efficiently  optimizing  the  orbital  coefficients ,{cul,  and  the  Cl 
coefficients  (c^l  simul taneously. ^3-39  -j^is  approach  exploits  an  idea  due  to 

/  35) 

Levy  of  using,  unitary  exponential  operators  to  perform  the  orbital 
rotation.  The  same  approach  can  be  used  within  the  coupled-cluster  framework 
which  will  have  the  advantages-,  (1)  that  all  results  are  size-extensive;  and 
(2),  that  the  configuration  expansion  is  not  just  a  few  determinants,  but 
consists  of  all  single,  double,  and  (via  the  nonlinear  coupled  cluster  scheme) 
the  quadruple  excitations. 


i 
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Consider  the  Hamiltonian  H,  related  to  the  usual  electrostatic 
hamiltonian  by 

H  =  U+HU  (1) 


where  U  is  a  similarity  transformation,  which  we  will  choose  to  be  unitary. 
The  uni tarity  c?r.  be  ensured  by  defining 

U  =  ein  =  e*  .  (2) 

For  n  a  hermitian  operator,  which  makes  k  skew  hermitian  (i.e.,  =  -k). 


where 


Lt  “1C,.  < 

n  =  e  He 


(3) 


K 


£  K  (x+x  -  X*X  ) 

r>s  rs  r  s  5  r 


(4) 


The  eigenvalues  of  any  hamiltonian  of  the  form  in  Eqs.  (1)  and  (3)  are 
unchanged  by  the  transformation,  so  this  transformed  H  may  be  treated 
just  as  the  ordinary  hamiltonian,  although  an  additional  degree  of  flexibility 
is  introduced  by  the  transformation. 

Employing  the  coupled-cluster  ansatz. 


fcc  =  e  IV 


T  =  T]  +  T2  +  ... 


T  =  1/nlrtjE***  xW...xv.X, 

n  ijk...  abc  kji 


(5) 


W 
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for  <p  a  single  determinant  reference  function  and  ^  intermediately 
normalized,  we  can  back-project  ^  onto  a  sufficient  set  of  single,  double, 
etc.,  excitations  to  give  the  CC  equations. 


a.  +  y  b..t.  +  y  c  • .  1 1 .  t,  +  y  +  +  +  _ « 

1  j  1J  J  j,k  1jk  J  k  j , k ,  1  dijkl  jW" 


(6) 


for  the  various  amplitudes  f t^ }  in  the  operator  T.  The  quantities  a., 

b..,  c... ,  etc.  are  combinations  of  integrals  relative  to  the  molecular 
i  J  l  j  k 

orbital  basis  set. 

If  we  limit  ourselves  to  T=T2,  then  we  obtain  the  coupled  cluster 

doubles  (CCD)  model,  which  terminates  Eq.  (6)  after  quadratic  terms,  and 

requires  that  we  determine  as  many  amplitudes  t.  =  tab  as  there  are  distinct 

double  excitations.  If  we  approximate  T=T^  +  T2,  then  T-|  will  occur  to  the 

2 

fourth-power,  T 2  quadratically,  and  the  coupling  terms  and  T2  will 

also  contribute.  In  this  case,  we  have  both  the  single  excitation  amplitudes, 
t®,  and  double-excitation  amplitudes  to  obtain.  The  energy  is  given  from  the 
Schrondinger  equation  by 

E  *  <*0M<tc>  ^oi^lV  »  (7) 


=  !  <ij||ab>  t*b.  +  l  <i|h,|a>ta 
i>J  i,a  1  1 


+  I  l  <ijl!aj>ta  +1 12  l  l  <ij  |  jab>  tat*? 
i,a  j  1  i,a  j,b  1  J 


H 


I 
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or 


l 

i>j 

a>b 


<1j| |ab>(t?j  +  i » j)  + 


l 

i,a 


<i 


!  f ! a  >  tj 


(8  ) 


where  the  usual  Fock  operator  expression  is  introduced. 

Eq.  (8)  provides  the  energy  given  the  amplitudes  for  double 

excitations,  t^,  and  single  excitations  t^  determined  by  Eq.  (6).  The 

molecular  integrals  and  amplitudes  in  Eq.  (8)  pertain  to  the  molecular 

orbitals  associated  with  the  transformed  hamiltonian,  H .  These  orbitals 

are  permitted  to  change  to  assume  a  more  optimum  form  by  exploiting  the 

transformation  in  Eq.  (3).  Since  the  coupled  cluster  results  are  invariant 

to  any  unitary  transformation  that  only  mixes  excited  orbitals  (i.e., 

a,b,c,d,etc. )  or  occupied  orbitals  (i.e.,  i,j,k,l,  etc.)  among  themselves, 

the  changes  due  to  the  transformation  arise  by  mixing  the  occupied  orbitals 

with  the  excited  orbitals.  The  transformation  matrix  {*  }  can  be 

rs 

determined  from  the  Hansdoff  expansion  by  imposing  a  stationary  condition 
on  the  orbitals. 


From  Eq.  (7)  and  Eq.  (3),  we  have 

E  =  <*0|e"KHelc|#cc>  .  (9) 

Invoking  the  Hausdorff  expansion  through  terms  second-order  in  k, 


(10) 


1 
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The  stationary  condition  is  accomplished  by  varying  E  with  respect  to  the 
operator  k.  Hence, 

5  E  =  0  =  <4>0 1  [H,  6<J  l’^cc>  +  1/2  <4>0I  [  [H,«k]k]  Kcc> 

(ID 

+  1/2  <$0!  [  [H,k]Sk]  kcc> 

Eq.  (11)  defines  the  matrix  elements,  <rs*  in  the  operator  <,  via 

T-— —  =  0  =  <4>  I  [H ,  X+X  -X:X  ]Urr> 

3*c  *o!  L  ’  r  s  s  rJ,vCC 

rs 

+  1/2  <*0|[[H,X*X$  -  X^Xr]<] l^cc>  (12) 

+  1/2  <*Q|  [[H,<c](xJXs  -  XsV)]|^CC> 

Since  <  appears  linearly  in  Eq.  (12),  by  working  out  the 

commutator  expressions,  an  equation  for  the  matrix  elements  {<  }  is 

obtained.  With  these  quantities  defined,  the  transformation  matrix  in 

Eq.  (1),  U  =  e<,  is  determined.  Hence,  Eq.  (12)  provides  k  for  stationary 

values  of  the  energy  (assuming  the  transition  state  form).  This  defines 

a  new  set  of  molecular  orbitals,  $  =  j  U  .  All  molecular  integrals  are 

u  L  v  vu 

transformed  to  this  new  set.  Then  the  coupled-cluster  equations,  Eq.  (6) 
can  be  solved  again  for  the  transformed  orbitals.  Successive  repetition 
of  this  procedure  should  give  a  coupled-cluster  solution  with  an  improved 
molecular  orbital  basis  set. 
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Since, 

*CC  =  eTlV  =  1  +  Tl+T2  +  1/2T1  +  1/2T2  +  T3  +  l*0>  (13) 

the  transition  state  formula  of  Eq.  (10)  only  permits  T-j  and  7^  to  occur 
linearly  in  the  determination  of  <,  since  the  other  terms  would  have  vanishing 
matrix  elements  with  |<j>  >.  This  seems  like  a  reasonable  approximation  for  most 
problems,  but  in  pathological  cases  the  quadratic  and  higher  terms  become 
increasingly  impor*  nt,  and  it  is  just  such  pathological  cases,  where  the 
orbital  optimization  is  most  useful. 

Alternatively,  variational  expressions  of  the  form 

E  =  ^^le  lcHe>c ] ^cq>/ <^cc I H4) 

can  be  employed  to  determine  equations  for  the  { icrs >  matrix  elements. 

This  has  the  advantage  that  a  rigorous  upper-bound  expression  is  , 
so  the  orbitals  should  be  "better",  and  that  higher  powers  of  T  operators 

would  contribute  to  this  determination.  But  it  is  still  necessary  to  truncate 
the  Hausdort'f  expansion.  Hence,  the  method  is  at  best  quasivariational . 

Furthermore,  it  is  not  clear  at  just  what  levels  appropriate  for  the 
trunction  of  T  or  <.  Additional  problems  pertain  to  T  being  a  non- 
hermitian  operator  while  <  has  to  be  chosen -to  be  hermitian  to  maintain 
the  unitary  transformation  property,  and  that  unlike  the  usual  coupled-cluster 
equations,  whose  commutation  using  the  normal  T,  i.e.,  e'THeT|$o>  terminates  at 
four-fold  commutations,  a  unitary  T  or  k  operator  requires  an  infinite 


commutator  series. 
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We  are  currently  investigating  the  relationships  among  the  several 
different  ways  of  imposing  variational  conditions  on  the  orbital  optimized 
coupled-cluster  equations. 
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ABSTRACT 
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A  goal  of  quantum  chemistry  in  biomedical  sciences  is  to  provide 
accurate  calculations  of  molecular  interaction  among  biochemical  molecules, 
drugs,  carcinogens,  etc.  In  this  effort,  there  is  a  natural  progression  from 
semi -empirical  quantum  chemistry,  to  ab  initio  self-consistent  field  theory, 
to  methods  that  properly  include  electron-correlation.  As  ab  initio  theories 
continue  to  develop,  many  more  problems  of  biomedical  interest  can  be 
addressed  by  accurate  correlated  methods.  The  intent  of  this  contribution  is 
to  discuss  many-body  approaches  to  the  correlation  problem,  i.  e. ,  many-body 
perturbation  theory  (MBPT )  and  coupled-cluster  methods  (CCM).  Unlike  most 
configuration  interaction  (Cl)  methods,  MBPT/CCM  offersa  number  of  important 
features  in  the  extension  to  larger  molecules.  These  include  the  proper 
dependence  of  the  correlated  calculation  on  the  size  of  the  molecule  (i.e., 
size-extensivity) .  This  has  significant  consequences  for  predictions  of 
ground  and  excited-state  properties.  These  features  will  be  illustrated  by 
applications  to  selected  molecules.  It  will  be  demonstrated  that  MBPT/CCM 
offers  a  natural  generalization  of  SCF  theory  that  is  formally  suitable  for 
applications  to  some  of  the  molecules  that  occur  in  biomedical  studies. 


I. 


INTRODUCTION 


In  the  applications  of  quantum  chemical  methods  to  problems  involving 
biochemical  molecules  and  their  interactions,  there  is  a  natural  progression 
from  empirical  or  semi -empirical  models  and  methods  to  ab  initio  self-consistent 
field  (SCF)  approaches,  and  eventually,  to  ab  initio  approaches  that  properly 
include  the  effects  of  electron  correlation.  The  purpose  of  this  contribu¬ 
tion  is  to  discuss  the  many-body  methods (i.e.,  many-body  perturbation 
theory^1'3}  [MBPT]  and  coupled  cluster  methods^4"7^  [CCM])  for  including 
electron  correlation.  The  emphasis  is  on  the  advantages  that  these  methods 
offer  over  the  more  traditional  Configuration  interaction  (Cl)  approaches  in 
large  molecule  applications. 

Semi-empirical  models  and  methods,  which  should  be  used  synergistical ly 
with  experiment,  are  most  properly  employed  to  investigate  trends  among  a 
series  of  similar  molecules.  Such  methods  can  be  used  for  rather  large 
molecules  relatively  inexpensively,  and  are  thus  finding  wide  use  in  bio¬ 
chemistry  and  particularly  in  drug  design.  on  the  other  hand,  in  principle 

ab  initio  methods  can  provide  hard,  quantitative  results  for  molecular  systems, 
which  can  be  potentially  used  to  complement  various  experimental  methods  by 
providing  answers  to  classes  of  problems  that  are  not  as  amenable  to  experiment. 

An  example  would  be  identifying  the  transition  state  and  activation  barrier 
in  a  reaction. 

In  practice,  however,  ab  initio  quantum  chemistry  suffers  from 
severe  limitations,  that  have  only  permitted  highly  accurate  results  to  be 
obtained  for  comparatively  small  molecules.  These  limitations  are  basically 
of  three  types:  (1)  Number  of  degrees  of  freedom  in  molecular  systems; 


(2)  limited  size  of  basis  set  that  can  be  used;  and  (3)  required  degree  of 
accuracy  of  the  ab  initio  approach. 

In  the  first  category,  the  problem  essentially  revolves  around 
the  Born-Oppenheimer  (or  fixed  nuclei)  approximation,  since  the  calculation 
of  the  electronic  structure  and  energy  must  be  repeated  for  each  choice  of 
coordinates  for  the  nuclei.  Limitations  (2)  and  (3)  pertain  to  each  of  these 
calculations  while  limitation  (1)  refers  to  the  number  of  times  the  calcula¬ 
tions  must  be  repeated.  For  example,  mapping  out  a  potential  energy  surface 
for  even  a  four-atom  system  with  3N-6=6  degrees  of  freedom,  and  computing 
10  points  for  each  degree  of  freedom,  would  amount  to  a  million  calculations. 

In  quantum  biochemistry,  fortunately,  one  is  not  often  interested  in  a  complete 
energy  surface,  but  usually  only  a  few  crucial  bond  lengths  and  angles  that 
need  to  be  optimized,  but  this  is  still  a  formidable  problem.  The  development 
of  SCF^”^  and  correlated  gradient  methods is  a  welcome  addition  to 
the  quantum  chemist's  repetoire,  but  even  these  techniques  are  only  applicable 
to  a  few  degrees  of  freedom. 

To  take  an  example  in  quantum  biochemistry,  consider  a  solvated 
molecule  where  it  is  recognized  that  the  solvation  characteristics  are 
partially  responsible  for  the  conformation  of  the  molecule  which  can  directly 
affect  a  highly  specific  interaction.  The  only  feasible  approach  to  such 
a  problem  at  present  is  the  determination  of  analytic  model  potentials  of  the 
Lennard-Jones ,  generalized  Morse,  and  other  types,  with  parameters  deter- 
mined  empirically,'  or  perhaps  from  highly  accurate  quantum  chemical 
calculations  of  the  component  pieces  of  the  larger  system. (28»29)  Then, 
these  potentials  can  be  used  to  handle  most  of  the  dynamical  movement  of  the 
molecule  and  solvent,  allowing  the  more  accurate  quantum  chemical  methods 
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augmented  by  gradient  techniques  to  focus  on  the  most  crucial  active  site 
interactions.  The  results  of  this  procedure,  however,  are  no  better  than 
the  accuracy  of  the  individual  calculations  which  are  subject  to  limitations 
(2)  and  (3). 

In  Figure  1  is  shown  a  schematic  drawing  that  illustrates  the 
dependence  of  an  ab  initio  quantum  chemical  prediction  on  basis  set  and 
calibre  of  method.  Just  improving  the  basis  set  or  method  is  not  enough, 
but  rather  a  systematic  improvement  in  both  is  required. 

Considering  the  basis  set  problem  first,  and  depending  upon  the 
property  of  interest,  it  is  a  matter  of  opinion  at  just  how  many  basis 
functions  are  required  to  obtain  good  SCF  results  for  molecules;  but,  certainly 
one  would  want  at  least  a  minimum  basis  set  of  one  Slater  orbital  (or  contracted 
Gaussian  orbital,  i.e.,  SZ  -  single  zeta)  for  each  electron  and  probably  two 
(DZ,  double  zeta)  or  more  (DZP,  double  zeta  +  polarization).  The  number  of 

4 

molecular  integrals  needed  to  do  an  SCF  calculation  rises  formally  as  n 

where  n  is  the  number  of  basis  functions  although  for  sufficiently  large 

2 

molecules  this  dependence  can  be  reduced  to  n  .  The  largest  SCF  calculations 
which  have  been  done  employ  no  more  than  ^300  functions.  This  imposes  a 
limit  of  at  most  300  electrons,  or  more  realistically  'vl 00  to  150  electrons 
explicitly  considered. 

The  problem  is  further  compounded  when  electron  correlation  is 
included.  Except  for  second-order  perturbation  theory  which  will  be  considered 
in  more  detail  in  Section  IV,  correlated  methods  have  a  dependence  on  the 
number  of  basis  functions  of  'wi^.  Again,  it  is  possible  to  reduce  this  by 


maybe  two-orders  of  magnitude  for  sufficiently  large  molecules,  but  it  is 
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evident  that  even  fewer  problems  can  be  studied  at  the  correlated  level  than 


at  the  SCF  level . 


(30,31) 


There  has  yet  to  be  a  really  good  idea  for  eliminating  the  basis 

set  problem  in  quantum  chemistry.  Completely  numerical  SCF  calculations 

(32) 

have  only  been  accomplished  for  a  few  diatomic  molecules,  and  nothing 
of  general  utility  has  yet  emerged.  At  the  cost  of  using  unrealistic 
potentials,  the  numerical  procedures  of  the  type  used  in  MS-X  have  had  some 

a 

(33) 

success.  From  the  viewpoint  of  basis  set  quantum  chemical  computations 

the  development  of  effective  potentials  for  the  chemically  inert  electrons 

(34-38)  (39) 

in  heavy  atom  molecules  is  very  useful.  '  Also,  Gaussian  lobe  functions 

chosen  to  represent  the  bonds  in  a  molecule  rather  than  located  on  the 

various  atomic  centers  have  reduced  the  number  of  basis  functions  while 

simplifying  the  calculations  of  the  integrals.  Various  integral 

(31 ) 

approximations'  '  and  other  clever  schemes  can  also  aid  in  making  the  calcula¬ 
tions  more  efficient,  but  the  basis  set  problem  remains  a  fundamental  limitation. 

The  third  limitation  above,  as  illustrated  in  Figure  1,  pertains  to 
the  degree  of  accuracy  of  method  that  is  required  for  the  property  of  interest, 
which  is  the  main  concern  of  this  contribution.  Generally,  SCF  theory  is 
considered  to  be  adequate  (+  10%)  for  molecular  conformations,  equilibrium 
molecular  structure,  and  first-order  properties,  i.e.,  properties  obtained 
as  an  expectation  value  over  the  SCF  density,  such  as  the  electrostatic 
potential  or  dipole  moments.  Oh  the  other  hand,  correlated  methods  are 
considered  absolutely  necessary  to  predict  electronic  and  photoelectronic 
spectra,  to  study  binding  energies  and  other  thermochemical  quantities  in 
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reactions  where  bond  breaking  is  occurring; and, most  second-  and  higher-order 
properties  like  polarizabilities,  sheilding  constants,  magnetic  succeptibilities, 
etc.  Since  many  questions  in  quantum  biochemistry  revolve  around  one  or 
another  of  the  properties  that  need  an  accurate  treatment  of  electron  correla¬ 
tion,  it  is  important  to  consider  the  characteristics  that  a  correlated  method 
should  have  if  it  is  to  be  applied  to  the  large  molecules  that  occur  in 
quantum  biochemistry. 

A  few  desirable  characteristics  for  such  a  correlated  approach 
are  that  the  method  should  be 

•  size-extensive*  (i.e.,  should  scale  properly  with  the  size  of  molecule) 

•  generally  applicable  to  a  wide  class  of  problems;  (i.e.,  avoid 
specific  formulations  or  choices  of  configurations.) 

•  efficient  and  cost-effective  (i.e.,  provide  large  correlations 
corrections  inexpensively) 

•  applicable  to  open-shells  and  excited  states; 

•  able  to  correctly  separate  a  molecule  into  its  fragments. 


Another  condition  that  one  might  expect  is  that  the  method  be  variational, 
giving  an  upper  bound  for  the  total  energy.  Lacking  a  coordinate  lower  bound, 
we  believe  this  is  an  unnecessary  restriction  since  the  quantities  of  interest 
in  quantum  chemistry  are  invariably  energy  differences  like  binding  energies, 
which  possess  no  rigorous  variational  properties  even  if  the  individual 
calculations  are  variational.  Furthermore,  except  for  a  full  Cl,  and  a  few 
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other  isolated  cases  (generalized  valence  bond,  GVB^4^,  e.g.)  a  variational 
requirement  is  not  consistent  with  the  size-extensive  condition  above,  which 
is  felt  to  be  much  more  important  to  satisfy  for  large  molecules. 

At  the  present  state-of-the-art  in  correlated  theory,  the  first 
three  conditions  are  easily  accomplished  with  MBPT.  Any  approach  based  upon 
the  linked-diagram  theorem  is  size-extensive.  A  large  class  of  problems  can 
be  studied  with  single  reference  MBPT/CCM  calculations  provided  that  RHF 
(of  a  UHF  open-shell  solution)  is  an  adequate  starting  point. 

For  the  cost-effective  property,  it  will  be  shown  that  second-order  perturba¬ 
tion  theory,  which  is  the  simplest  MBPT  approximation,  typically  accounts 
for  ^90  percent  of  the  correlation  energy  in  a  basis  set  and  significantly 
improves  the  SCF  predictions  of  dissociation  energies  and  molecular  geometries. 
Since  this  requires  only  marginally  more  effort  than  an  SCF  calculation,  since 
it  is  size-extensive, and  has  rather  general  utility,  it  is  a  very  attractive 
lowest-order  approximation. 

The  fourth  requirement  can  be  handled  with  many-body  approaches  such 
as  equation-of-motion  techniques ,(42-44)  or  wl-th  Cl,  and  the  fifth 
is  currently  most  easily  achieved  using  Cl  methods.  In  the  last  case,  the 
MBPT/CCM  theory  exists  for  this  problem,  but  has  not  yet  been  implemented 
in  a  general  purpose  program. (^-47)  jn  many  cases>  a  (jHF  reference  func¬ 
tion  will  permit  correct  separation,  but  the  path  toward  the  separated  limit 

(91 

is  not  always  accurate/  ' 
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In  Section  II,  the  size  extensive  property  of  MBPT/CCM  will  be 
discussed  in  some  detail  since  this  is  an  extremely  important  condition 
for  potential  applications  of  correlated  methods  to  large  molecules. 

Section  III  will  present  a  brief  discussion  of  the  ideas  in  many-body  theory 
that  are  important  for  large  molecules,  while  the  final  section  will  focus  on 
some  applications  to  benzene  to  demonstrate  the  nature  of  correlation  effects 
due  to  higher  excitations  in  this  prototype  system.  In  this  section  emphasis 
will  also  be  placed  on  the  accuracy  of  the  simplest  approximation,  second-order 
perturbation  theory,  which  typically  provides  a  very  large  part  of  the  electron 
correlation  effect  as  an  inexpensive  by-product  of  the  SCF  calculation. 


•H. 
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II. _ SIZE-EXTENSIVITY  IN  MOLECULAR  CALCULATIONS 

Probably  the  best  way  to  illustrate  the  importance  of  quantum 

mechanical  methods  that  scale  properly  with  the  size  of  a  molecule  is  to 

consider  the  model  problem  of  a  lattice  of  separated  electron  pair  bonds, 

such  as  H2  molecules,  since  this  serves  as  a  first  approximation  to  any  large 

(48  49  111 

molecule.  This  problem  has  been  worked  out  by  several  investigators/  *  ’  ' 

but  it  is  pertinent  enough  to  the  discussion  that  it  is  worth  presenting  a 
version  here. 

Assume  the  ^  molecules  are  sufficiently  far  apart  or  separated 
by  barriers  so  that  they  can  be  considered  to  be  noninteracting.  For  sim¬ 
plicity,  we  will  further  assume  that  the  component  set  of  molecular  orbitals 
for  each  molecule  are  natural  orbitals  so  that  single  excitations  in  the 
H2  wavefunction  need  not  be  considered.  Then  the  intermediately  normalized 
wavefunction  for  each  molecule  i,  may  be  written. 


Ai)  =  Ai)  +  xM(i) 


0) 


where  ( i )  is  the  first  natural  determinant  (close  to  the  Hartree-Fock 

solution)  and  X  (i)  is  a  sum  of  doubly-excited  determinants  including  their 
appropriate  coefficients.  The  norm  of  the  function  in  Eq.  (1)  is  1  +  S,  where 


S(i) 


M,  .  N, 

=  <x  (i) 


(2) 


<*q(i) |xM(i )>  =  0 


(3) 
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The  wavefunction  for  the  lattice  is 

4.  =  TT  (*  (i)  +  x  (i ) )  •  (4) 

i=l 

Antisymmetry  is  disregarded  since  the  molecules  are  noninteracting. 

With  H*-  =  X  H ( i ) ,  the  energy  of  the  lattice 


EL  = 


<<fLlHL|$L> 
<$  |  $  > 


=  NE, 


(5) 


where  is  the  energy  of  the  H2  molecule.  With  b  =  <*^(1) | H ( i ) |x^(i )>. 
which  is  essentially  the  correlation  energy  of  the  molecule, 


+  6 


(6) 


A  method  is  said  to  be  "size-extensive"  if  the  total  energy  calculated  by 
the  method  is  appropriately  linear  in  N,  as  in  Eq.  (5). 

Notice  that  the  product  wavefunction  in  Eq.  (5)  includes  quadratic 
and  higher  product  terms  like  xn(i)x n(j)»  which  correspond  to  simultaneous 
double-excitations  on  different  centers,  but  are  quadruple  and  higher- 
excitations  in  a  super-molecule  Cl  description.  Since  these  terms  arise 
from  disjoint  double  excitations,  they  are  fundamentally  simple,  but  the 
standard  Cl  framework  is  not  able  to  exploit  this  inherent  simplicity.  This 
causes  an  innate  error  in  truncated  Cl  that  becomes  most  important  for  larger 
molecules. 
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To  investigate  this,  we  can  consider  a  reference  wavefunction  for  the 
lattice  of  the  form 


with  energy 


N 

EL  =  l  EM( i )  =  NEM 
0  -j  — i  o  o 


(7) 


(8) 


Using  this  reference  function,  the  double-excitation  Cl  ( D-C I )  wavefunction 
for  the  lattice  is  constructed  as 

L  L  ^  L 

Vci  =  To  +k^  V  (k) 

=  *0(1)  *0(2)...*0(N)  +  ^  ck  ^xM(k)/$o(k)  (9) 

L  N  L 

■  +  c  l  r(k)  . 

0  k=l 


The  weighting  coefficient,  c,  is  the  same  for  each  molecule  in 
a  noninteracting  lattice.  Using  the  expressions, 


<xM(i)|xM(j-)  =  s^S 
<xM(i)|HL|xM(J)>  =  «ijtEH2s-e] 

ae  -  el-nem 
0 


(10) 
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it  follows  that 

<^|HL|yL(k)>  =  e  (11) 

<4'L(K)|HL|yLU)>  =  «k£CNSE^  +  (S-l)B]  .  (12) 

From  these  matrix  elements,  the  D-CI  secular  equation  becomes 

aE  -  Ncb  =  0  (13a) 


{B  -  [SaE  +  (S-l )b]}c=  0 


Solving  Eq.  (13)  simultaneously  for  aE  and  c. 


aEd-ci  =  8 


-(1-S)±  /fl-s)2  +  4SN ) 
2S 


Neffe  +  ism 


(13b) 


(14) 


The  positive  sign  is  required  since  aE<o  and  s<o.  Since  the  correct 
aE  *  Nb,  D-CI  is  not  size-extensive. 

With  the  aid  of  a  value  for  S  in  Eqs.  (2)  and  (14)  it  is  possible  to 
get  some  feel  for  the  size  of  these  effects.  From  a  natural  orbital  study  by 
Davidson  and  Jones^50^  of  the  50  term  Kolos-Roothaan^51 ^  wavefunction  for  Hj 
at  R  =  1.40,  S  for  ^  is  0.0181.  Some  representative  values  are  shown  in 
Table  I,  along  with  values  for  a  lattice  of  He  atoms  for  comparison  (SHe=0.0083) . 
It  is  apparent  that  the  error  in  the  correlation  energy  as  determined  by  D-CI 
can  be  significant  even  for  modest  numbers  of  electrons.  It  is  also  apparent 
that  the  errors  are  greater  for  typical  covalent  bonds  than  inner-shell  electron 


pairs  as  in  He  atoms.  In  fact,  we  will  find  that  Table  I  can  provide 
a  rather  accurate  estimate  of  the  effects  of  higher  excitations  simply 
by  counting  the  number  of  electrons  in  covalent  bonds  and  inner-shell 
electron  pairs. 

M  M 

Since  the  product  terms  x  (i )x  (j )  correspond  to  quadruple 

excitations  in  a  super-molecule  Cl,  while  triple  products  are  hextuple 

excitations,  etc.,  the  size-extensive  property  of  MBPT/CCM,  that  is  a 

(1  21 

consequence  of  the  linked-diagram  theorem,'  ’  '  is  essentially  a  result  of 

a  more  proper  treatment  of  quadruple  and  higher  excitations  than  in  Cl. 

Hence,  a  statement  that  size-extensivity  is  important  in  correlated 

calculations,  is  equivalent  to  the  statement  that  quadruple  and  higher- 

excitations  are  important.  Since  the  number  of  configurations  in  Cl  are 

proportional  to  the  number  of  basis  functions  raised  to  the  level  of 

excitation  included,  the  number  of  quadruple  excitations  generated  from  100  basis 

4  8 

functions  would  require  'v(lOO)  or  10  configurations.  Hence,  better 
computational  methods  for  including  effects  of  higher-excitations  in 
correlated  calculations  are  extremely  important.  Many-body  methods  tend  to 
take  the  intelligent  viewpoint  that  removing  the  erroneous  terms  (i.e.,  unlinked 
diagrams)  in  D-CI  is  preferable  to  including  higher-excitations.  In  practice, 
this  viewpoint  leads  to  computationally  more  tractable  equations  that  are 
closer  to  those  in  D-CI. (7) 

Since  the  correct  density  is  in  error  by  the  neglect  of  the  product 

(4°1 

terms  in  the  wavefunctions ,  further  analysis  discussed  elsewhere,  '' 
demonstrates  that  the  density  matrix  obtained  from  a  truncated  Cl  reduces  to 
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just  the  density  computed  from  the  reference  function.  If  the  latter  is 
an  SCF  function  we  have 


Lim  TCI  .  SCF 
N-*»  p  '  p 


(15) 


Similarly,  for  an  excitation  energy. 


Lim 

N-*» 


)  =  E 


(49) 

SCF 

f 


escf 

l 


(16 ) 


Hence,  size-extensi vity  affects  more  than  the  total  energy. 

One  additional  consequence  worth  mentioning  is  that  in  a  typical 
reaction, 

A  +  8  +  C  +  D  (17) 

the  heat  of  the  reaction,  AHrxn  *  AH^C)  +  aH^-(D)  -  AHf(A)  -  AHf(B). 

However,  if  these  individual  quantities  are  determined  by  a  truncated  Cl 
this  simple  addition  is  not  entirely  justified,  since  the  truncated  Cl 
ignores  the  simultaneous  excitations  that  prohibit  AHf(C+D)  at  from 

being  AHf(C)  +  AHf(D).  In  practice,  this  frequently  requires  that  one 
compute  the  super-molecules  C  +  D  and  A  +  B  in  Cl  to  make  the  energy 
difference  most  accurate.  This  should  be  contrasted  to  predictions  made  with 
a  size-extensive  method  where  a  table  of  results  for  species  obtained  at  a 
given  level  of  approximation  may  be  added  and  subtracted  just  like  the 
experimental  values. 


» 
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III. _ SYNOPSIS  OF  MANY-BODY  THEORY 

The  theory  of  MBPT/CCM  has  been  discussed  in  detail  in  several 
places. in  particular  reference  [7]  provides  a  fairly  detailed  mathe¬ 
matical  description  from  the  viewpoint  taken  in  this  article.  The  theory 
as  originally  developed,  uses  second-quantization  and  diagram  techniques, 
which  are  unfamiliar  to  many-quantum  chemists,  and  this  tends  to  camouflage 
the  important  concepts  that  emerge  from  the  many-body  approach.  Instead  of 
presenting  any  detailed  mathematical  development  here,  we  will  sketch  the 
basis  for  the  two  significant  concepts  that  emerge  from  MBPT/CCM,  namely 
the  linked-diagram  theorem,  which  guarantees  size-extensivity,  and  the  cluster 
decomposition  of  Cl  excitations  into  separate,  more  physically  satisfying 
pieces  that  lead  to  tractable  equations  for  including  the  effects  of  higher 
excitations.  Consult  reference  [7]  for  detailed  equations,  and  the  original 
references'”6^  for  the  complete  formal  development.  For  simplicity,  in  the 
following  we  will  limit  ourselves  to  a  single  reference  function  such  as  an 
unrestricted  Hartree-Fock  (SCF)  solution.  Various  versions  of  the  multi- 

(45-47) 

reference  function  theory  are  available. 

It  is  well-known  that  one  way  to  solve  a  Cl  eigenvalue  equation 

(52) 

is  with  perturbation  theory.  '  Using  the  Rayleigh-Schrodinger  form,  we 
can  separate  the  hamiltonian  H  =  HQ  +  V,  where  H0  is  the  sum  of  the  SCF 
one  electron  hamiltonians  and  V  is  the  two-electron  part  minus  the  SCF  effective 
one  particle  hamiltonian,  then  we  have 


*  m 
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H  <f> 
oo 


=  E 


oTo 


II  nflbc  ■  •  _  p  3  bC  «  •  ndbc  • 

H0  °ijk”  ■  Ei jk‘  •  Dijk- 


(18) 


dbc 

for  4>  the  SCF  solution,  and  the  various  determinants  that  can  be 

o  *  jk  •  •  • 

formed  by  replacing  occupied  SCF  orbitals  with  excited  SCF  orbitals.  The 
Cl  eigenvalue  through  fourth-order,  then  becomes 


E  •  E„*  <<(JK|*0>  *  <t0|VR0Vk0>  +  <*0|VR0(V-<V.)R0V|V  ()9) 

*  <*0|VR0(V-<V>)  R0(V-<V>)  R0V|*0.  -  <*0|VR0V|V<*0|VR|V|V  . 

The  resolvent,  RQ,  has  the  expression 

R0  =  !(L><i!lEo"H0,->"1<^  {20) 

where  |h>  is  composed  of  all  the  Cl  excitation  .  Even  though  |h> 

is  formally  complete, Slater's  rules  for  matrix  elements  chooses  from  all 
possibilities  only  the  few  that  have  nonvanishing  contributions- 

Subject  to  an  SCF  reference  function,  only  double-excitations  can 
mix  across  V  with  <ji0,  so  the  second-  and  third-  order  terms  in  Eq.  (19) 
involve  only  double-excitations.  The  first  of  the  two  terms  in  fourth-order, 
however,  can  mix  single,  double,  triple,  and  quadruple  excitations  at  the 
middle  Rq,  although  the  second  fourth-order  term  (i.e.,  the  renormalization 
term)  still  has  no  contributions  except  from  double  excitations. 

From  the  model  problem  of  separated  H2  molecules  presented  in 
the  previous  section,  it  is  easy  to  check  whether  each  of  the  terms  in 
Eq.  (19)  is  size-extensive.  Considering  the  second-order  energy  of  the  H2 
lattice 
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E2  =  <*0lVRoVl*0>  (21) 

as  an  example,  we  have  from  Slater's  rules 
i  nocc  nexc  , 

E2  *  1/4  .1.  I  l-iJl  |ab>|Z/(ei  +  e.  -  £fl  -  )  .  (22) 

i ,  j  a ,  b  J 

The  notation  <ij||ab>  =  <ij|ab>  -  <ij|ba>  =  (ta|bj)  -  (ib|aj),  and  {e^and  {e a } 
are  the  SCF  orbital  energies  for  the  occupied' and  excited  orbitals,  respectively. 
Using  a  little  algebra, 

^2  =  2  J  l  [(ia[jb)2  -  ( i a | j b )  (ib|ja)>(5i  +  £j  -  £a  -  ^  •  (23) 
i <j  a<b 

Since  this  expression  is  invariant  to  any  unitary  transformation  among 
degenerate  orbitals,  we  may  choose  the  orbitals  to  be  localized  on  the 

molecules  in  the  lattice  to  make  the  argument  most  transparent.  In  this 
case  the  only  nonvanishing  integrals  have  the  charge  distribution  (ia),  (jb), 

(ib)  or  (jc)  on  the  same  H2  molecule,  otherwise  the  terms  would  vanish. 

Hence,  it  follows  that 

Eg  =  NE^  (24) 

and  second-order  perturbation  theory  is  size-extensive.  It  can  be  similarly 


shown  that  this  is  also  true  for  E3. 
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Now  consider  fourth-order  .  The  renormalization  term  is  composed 
of  an  E^  term  and  a  similar  term  a  =  -<4)o  |  VR^ V  1 4>0>  which  differs  from  E2  only 
by  having  the  denominator  squared.  Since 

=  NE2M,  aL  =  ,  (25) 

2 

The  product  of  the  two  has  an  N  dependence,  which  is  erroneous.  If  E^  is 

to  be  size-extensive,  the  first  term  in  E^,  must  also  have  an  equal  and 

2 

opposite  N  dependence  to  cancel  out  these  uncharacteristic  terms.  The  single-, 
double-,  and  triple-excitation  contributions  to  the  first  part  of  E^can 
be  shown  to  be  size-extensive.  Hence,  to  resolve  the  problem,  it  is  necessary 
to  consider  the  quadruple  excitation  contributions.  Following  a  great  deal  of 
algebra!7  ^  the  quadruple  excitation  part,  E^,  may  be  written  in  the  form, 

=  E2A  +  Q  (26) 

where  Q  is  properly  size-extensive.  Hence,  E2a  cancels  the  renormalization 

2 

term  and  with  it,  the  erroneous  N  dependence.  This  is  the  substance  of  the 

linked-diagram  theorem.  The  algebraic  analysis  that  leads  to  Eq.  (26), 

represents  Q  as  linked  diagrams,  while  E2a  corresponds  to  unlinked  diagrams. 

A  similar  analysis  will  apply  in  all  higher  orders  which  is  the  linked- 
(  2  ) 

diagram  theorem.  This  provides  the  expression  for  the  energy, 

E  =  E„+  i  ^IVtVV'ViVL  |27> 

p=o 

where  L  limits  the  terms  to  only  1 inked-diagrams. 
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It  should  be  evident,  that  if  quadruple  excitations  had  not 

2 

been  included  in  E^,  then  the  term  with  the  erroneous  N  dependence  would 
remain.  This  is  exactly  what  happens  when  a  truncated  Cl  calculation  is  made. 
Limiting  the  configurations  to  single-  and  double-excitations,  for  example,  will 
necessarily  retain  these  erroneous  terms  destroying  the  size-extensivity  of 
the  method.  If  quadruple  excitations  were  to  be  included  in  the  Cl,  the  result 
would  be  size-extensive  through  fifth-order,  but  would  fail  in  six  and  higher- 
orders  due  to  hextuple  excitations.  On  the  other  hand,  any  approximation  to 
the  linked  diagram  theorem,  Eq.  (27),  is  size-extensive.  This  means  that  even 
second-order  perturbation  theory  can  be  much  better  than  very  good  Cl  calculations 
for  sufficiently  large  molecules. 

For  small  molecules,  multi -reference  Cl  techniques,  that  incorporate 
at  least  the  most  important  quadruple  excitations  as  double-excitations  from 
a  double-excitation  reference  space,  will  be  size  extensive  for  most  practical 
purposes.  GVB  calculations  are  size-extensive,  but  GVB-CI  will  be  only  approxi¬ 
mately  size-extensive  unless  all  excitations  into  the  GVB  orbitals  are  included. 
Since  GVB  provides  a  better  choice  of  orbitals  than  SCF,  and  since  one  includes 
higher-level  excitations  than  is  normal  in  SCF  based  Cl  approaches,  GVB-CI 
will  usually  be  closer  to  size-extensive  than  other  Cl  methods.  An  added 
advantage  is  that  within  the  GVB  method  it  is  often  possible  to  ensure 
correct  separation. 

The  other  important  idea  developed  in  many-body  theory  is  the 
cluster  expansion  of  the  wavefunction.  The  basic  concept  is  that  the  correct 
wavefunction  may  be  written  as  eTj<}>0>,  where  T  is  an  operator.  This  form  of 
the  wavefunction  ensures  the  linked-diagram,  size-extensive  basis  of  the  theory. 


A 
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Then  T  has  the  form 
T  =  T 


1  +  T2  +  T3  + 


(28) 


where  T-j,  Tj,  ...  are  one-body,  two-body,  ...  cluster  operators.  The  T^ 
operator  generates  double-excitations  with  amplitudes  to  be  determined  by 
the  coupled-cluster  equations,^  but  the  exponential  form 


eT  =  T  +  1/2T2  +  1/3! T3  +  ...  (29) 

causes  some  very  different  things  to  happen  than  in  the  Cl  approach.  To  see 
this,  consider  the  Cl  operator  for  quadruple  excitations,  C4<  By  equating 
the  Cl  and  coupled-cluster  expressions  for  the  quadruple  excitations,  we  have 


C4  =  T4  +  1/2T2  +  1/4  T.j4  +  1/2  T2  Tg  +  T^  .  (30  ) 

Physically,  what  does  this  mean?  Roughly,  T4  represents  an  interaction 
? 

among  four  electrons  while  T£  represents  two  simultaneous  interactions  of 

two  electrons.  A  transformation  to  Brueckner  orbitals  makes  T-j  vanish,  while 

T-j  is  usually  small  even  for  SCF  orbitals  so  the  final  three  terms  are 

negl igible  most  of  the  time.  Since  the  normal  electrostatic  hamiltonian  has 

no  more  than  two-electron  operators,  simultaneous  two-electron  interactions 

would  seem  to  be  much  more  frequent  in  molecules,  than  "true"  four  electron 

interactions.  From  another  viewpoint,  the  NF^  lattice  problem  emphasizes 

the  neglect  of  simultaneous  double-excitations  on  different  Hj  molecules, 

2 

which  is  exactly  what  T£  offers.  Thirdly,  from  perturbation  theory,  it  may 

2  (7) 

be  shown  that  all  the  fourth-order  quadruple  excitation  terms  arise  from  Tg,  ' 


i 


s-  < 


with  only  contributing  in  fifth-  and  higher-orders.  Consequently,  it 

was  suggested  by  Sinanoglu^1^  that  =  l/ZTg  is  a  very  good  approximation. 

Using  this  ansatz,  we  have  the  coupled-cluster  doubles  (CCD)  approximation 

T? 

for  the  wavefunction,  e  [$  >  •  This  leads  to  a  set  of  nonlinear  equations 
for  the  amplitudes  but  there  are  only  as  many  of  these  amplitudes  as  in 


a  D-CI 


(5,7) 


This  provides  the  benefit  that  we  have  a  size-extensive  method;  it 


is  infinite-order  although  restricted  to  T^,  and  we  have  no  more  amplitudes 
that  in  a  D-CI  calculation  even  though  the  effect  of  quadruple  excitations 
are  included.  Since  CCD  reduces  in  fourth-order  perturbation  theory  to  double- 
and  quadruple-excitation  diagrams,  it  is  straightforward  to  solve  the  CCD 
equations  as  successive  iterations  of  a  fourth-order  MBPT  calculation  (7). 
Hence,  couple-cluster  methods  may  be  viewed  as  complementary  to  MBPT  when 
higher-order  corrections  are  needed  as  can  become  important  in  pathological 
cases 


*1 
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IV. _ ILLUSTRATION  OF  MBPT/CCM  RESULTS  FOR  SOME  SMALL  MOLECULES 

The  simplest  approximation  to  the  correlation  energy  in  MBPT 

(assuming  an  SCF  reference  function  for  simplicity)  is  given  by  the  second- 

order  perturbation  theory  expression  of  Eq.  (22). 

Since  the  molecular  orbitals  i,  j,  ...a,  b,...  and  their  orbital 

energies  e- ,  e.,...e  ,  eK,  are  obtained  from  an  SCF  calculation,  all  the 

information  is  available  that  is  needed  for  a  correlated  calculation  except 

for  a  partial  integral  transformation,  if  used.  The  SCF  calculation  generates 

a  set  of  two-electron  integrals  relative  to  atomic  (i.e.,  primitive)  basis 

functions,  and  in  the  general  case  an  integral  transformation  is  required 

to  obtain  the  integrals  relative  to  the  molecular  orbitals,  i.e.,  <ab||cd>, 

5 

which  depends  on  the  number  of  basis  functions  as  n  ;  or,  alternatively,  a 

direct  calculation  of  E^,  E,,  in  terms  of  the  integrals  relative  to  atomic 

orbitals  (probably  orthogonal ized)  is  required.  In  the  case  of  E£,  however, 

only  a  very  small  number  of  integrals  are  required,  since  each  integral 

involves  only  two  occupied  and  two  unoccupied  orbitals.  Consequently,  E2 

2  2  4 

requires  no  more  than  n  n  <n  operations,  or  less  than  in  the  SCF  calculation 

OCt  6XC 

itself.  In  a  sufficiently  large  molecule  where  the  primitive  integrals  ( ctS ! y6 ) 

vanish  unless  a  and  s  are  in  the  same  neighborhood  as  are  y  and  5,  and 

unless  the  charge  distributions  (ag)  and  (y&)  are  not  too  far  apart,  the  SCF 

2 

calculations  goes  as  vi  .  In  this  case,  evaluating  E2  directly  in  terms  of 
(asjyi)  will  permit  a  similar  simplification,  hence  E2  can  always  be 
evaluated  as  a  by-product  of  large  SCF  calculations  at  negligible  additional 
expense. 

This  approximation  is  certainly  recommended  by  convenience,  but 
how  reliable  an  approximation  is  it  for  the  correlation  energy?  In  Table  II 
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are  shown  the  fractions  of  the  correlation  energy  within  a  basis  set  given  by 

E2,  E^,  and  the  fourth-order  contributions  just  from  double-  and  quadruple- 

(7  9) 

excitations  diagrams  *  for  a  variety  of  molecules.  Using  an  SCF  starting 

point,  E2  and  E^  are  solely  determined  by  double-excitations,  with  single-, 

double-,  triple-,  and  quadruple-excitations  appearing  in  fourth-order,  but 

in  the  interest  of  also  comparing  the  higher-order  corrections  obtained  by 

the  CCD  (coupled-cluster  doubles)  approximation,  which  includes  only  double- 

o 

excitations  and  the  disjoint  (i.e.,  T2)  quadruple-excitations  to  all  orders,  the 
single-  and  triple-excitations  contributions  are  omitted  from  Table  II. 

It  is  apparent  from  the  table  that  the  simple  second-order  approxima¬ 
tion  accounts  for  the  vast  majority  of  the  correlation  energy  obtainable 
within  the  basis  set.  A  few  general izations  about  the  results  may  be  made. 

In  multiply  bonded  systems  such  as  N2>  CO,  and  C02>  E2  tends  to  slightly 
overestimate  the  net  correlation  energy  in  the  basis  set,  while  for 
saturated  systems  like  H20,  CH^,  etc.,  it  is  more  likely  to  underestimate  the 
effect.  HCN  and  benzene  are  intermediate.  In  a  case  where  near  degeneracy 
plays  a  role  such  as  BH^,  convergence  of  the  perturbation  theory  is  much 
slower  making  E2  a  poorer  approximation.  No  particular  differences  are 
observed  for  open-shell  molecules  when  UHF-SCF  instead  of  an  RHF-SCF  solution 
is  used  as  the  unperturbed  solution.  On  the  average,  it  is  clear  that  E2 
accounts  for  %90  percent  of  the  correlation  energy  obtainable  within  the 

basis  set.  Since  these  basis  sets  are  good  enough  that  they  account  for  ^70 

(9) 

percent  of  the  "experimental"  valence  shell  correlation  energy;  this  means 

E2  gives  %60  percent  of  the  experimental  valence  correlation  energy.  It  is 

also  clear  that  DQ-MBPT (4 )  is  usually  very  close  to  the  infinite-order 
(9) 

CCD  model:  This  is  a  common  occurrence  except  for  cases  where  near  degeneracy 


is  a  problem. 
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In  Table  III  are  shown  some  thermochemical  results  obtained  from  E2 
compared  to  higher-order  correlation  approximations.  Although  E2  predictions 
are  not  as  good  as  the  better  approximations,  they  are  clearly  superior 
to  the  SCF  predictions,  again,  providing  most  of  the  observed  correlation 
corrections. 

A  similar  result  can  be  obtained  for  second-order  predictions  of 
molecular  structure,  where  on  the  average  -v50  percent  of  the  error  in  the 
SCF  geometries  is  removed. 

To  obtain  the  exceptional  accuracy  reflected  in  Table  III  and 
S3  541 

reported  elsewhere^  ’  '  for  various  properties  of  small  molecules,  it  is 

necessary  to  go  beyond  second-order,  but  for  large  molecules,  the  simplicity 
and  comparatively  high  accuracy  of  this  approximation  demands  that  it  be 
used  to  augment  any  large-scale  SCF  calculation  of  biochemical  interactions. 

Benzene  serves  as  a  prototype  of  many  of  the  large,  conjugated 
molecules  that  occur  in  biochemistry.  As  such,  it  is  appropriate  to 
analyze  the  higher  order  MBPT/CCM  description  of  electron  correlation  in 
benzene  to  develop  some  feeling  particularly  for  the  effect  of  quadruple 
excitations. 

The  basis  set  is  a  standard  Dunning  double  zeta  contraction  of 
Huuzinaga’s  9s5p  primitive  Gaussian  basis  for  carbon  and  the  two  H  Is  orbitals 
corresponding  to  a  Slater  exponents  of  1.2,  giving  72  CGTO.  The  SCR  energy 
of  -230.6369  differs  by  0.113  a.u.  from  the  SCF  results  for  a  DZP  basis  and  0.18 
a.u.  from  the  estimated  SCF  limit. ^8)  c(15^)  electrons  are  kept  frozen 

at  the  SCF  level  so  the  correlated  calculations  only  pertain  to  the  valence 
correlation  energy.  Polarization  functions  are  usually  found  to  be  more 
important  for  correlation  effects  than  in  the  SCF  calculation  itself,  so  the 
current  DZ  predictions  should  underestimate  the  magnitude  of 
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the  valence  correlation  energy.  Even  so,  it  is  apparent  from  Table  IV,  that 
quadruple  excitations  amount  to  ^20  percent  of  the  predicted  correlation  energy. 

In  an  attempt  to  study  the  origin  of  the  quadruple  excitation  effects, 
the  occupied  and  excited  pi-orbitals  were  removed  and  the  calculation  repeated 
to  give  a  value  for  just  the  sigma  framework  excited  solely  into  unoccupied 
sigma  orbitals.  The  same  procedure  was  carried  out  for  the  pi -electrons. 

These  results  are  reported  in  the  second  and  third  columns  of  Table  IV. 

The  sigma  framework  accounts  for  over  half  of  the  net  quadruple 
excitation  effect,  while  the  correlation  effects  of  the  delocalized  pi-electrons 
are  relatively  independent  of  the  quadruple  excitations.  Since  the  former 
involves  12  roughly  independent  covalent  bonds,  from  Table  I  and  Eq.  (14), 
one  would  estimate  an  effect  of  ^14  percent  in  reasonable  agreement  with  the 
calculated  12  percent.  The  effect  of  the  quadruple  excitations  on  the  pi- 
electron  bonds  is  much  smaller,  but  this  is  primarily  due  to  the  fact  that  only 
three  bonds  are  possible.  If  the  appropriate  S  for  the  pi -structure  were  as 
small  as  in  He,  the  estimated  effect  of  quadruple  excitations  would  be  1.6 
percent.  The  remaining  correlation  effects  come  from  the  sigma-pi  interactions. 
It  is  interesting  that  excitation  of  sigma  electrons  into  pi-excited  orbitals 
and  vice-versa  results  in  %13  percent  of  the  correction  energy. 

The  DZ  basis  used  here  is  capable  of  providing  only  about  49  percent 

of  the  experimental  valence  correlation  energy.  (58,59)  p0q arl- zati on  functions 

(9) 

would  improve  this  result  by  about  20  percent.  '  Since  the  quadruple  excitations 
are  also  responsible  for  more  than  20  percent  of  the  correlation  energy,  the 


to 


size  of  error  encountered  in  SD-CI  is  as  severe  as  excluding  polarization 
functions  from  the  basis  set.  Since  the  effect  of  quadruple  and  higher  Cl 
excitations  will  inevitably  increase  as  larger  molecules  than  benzene  are 
studied,  a  correlation  method  that  does  not  account  for  such  effects  has 
little  hope  of  contributing  to  the  investigation  of  biochemical  interactions. 
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In  every  example  but  CgHg,  the  basis  set  of  at  least  of  double  zeta  -plus  polarization 
quality.  In  a  double  zeta  basis  is  used. 
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TABLE  III.  COMPARISON  OF  THERMOCHEMISTRY  RESULTS 

OBTAINED  BY  SCF  AND  MBPT  WITH  EXPERIMENT 
[All  basis  sets  are  at  least  DZP  quality.] 


Reaction 


2  BH3-B2Hga 

SDQ-MBPT(4) 

BH3  +  C0-H3BC0a 

D-MBPT (4) 

BH3  +  NH3-H3BNH3a 

D-MBPT(4) 

HNC  -v  HCNb 

SDQ-MBPT (4) 

HNC^[HNC]b 

SDQ-MBPT(4) 

BNC  -  BCNb 

SDQ-MBPT (4) 

Li NC  -»-LiCNb 

SDQ-MBPT (4) 

CH3NC  -*CH3CNC 

SDQ-MBPT (4) 

CH3NC-[CH3"]C 

SDQ-MBPT(4) 

H  +  CO  -»■  HC0d 

CCD 

HCO  -*•  [HCO]d 

CCD 

h2co  -*-h2  +  coe 

CCD 

H,C0  -*-H  +  HCOe 

CCD 

-aE  (kcal/mole 


MBPT/CCD  Experiment 


36.6  +  2 


20.4  +  2 


(10.3  +  l)9 


23.7  +  14" 
-38. 41 

15.7  +  1.5 


-1.9* 

86.0  +  1.0 


>— n  rj 
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a  Reference  [53]. 

^  Reference  [55].  Square  bracket  indicates  a  transition  state. 

This  result  includes  a  4  kcal/mole  zero  point  correction 
for  the  transition  state. 

Reference  [54].  Square  bracket  indicates  a  transition  state. 

This  result  includes  a  4.8  kcal/mole  zero  point  correlation, 
for  the  transition  state. 

^  Reference  [56]. 

e  Reference  [57]. 

f  T.  P.  Fehlner  and  G.  W.  Mappes,  J.  Phys.  Chem.  73,  873  (1969). 

9  L.  Maki,  unpublished  results.  Reference  [55]  concludes  that  this  experi¬ 
mental  value  is  in  error.  The  result  should  be  15  +  2  kcal/mole. 

h  M.  H.  Baghal-Vayjovec,  J.  L.  Collister,  and  H.  0.  Pritchard,  Can.  J. 

Chem.  55,  2634  (1977). 

1  F.  W.  Schneider  and  B.  S.  Rabinovitch.  J.  Am.  Chem.  Soc.  65,  1794  (1969). 

J  P.  Warneck,  Z  Naturforsch  A26,  2047  (1971). 

K.  Yamada,  T.  Nagakuru,  K..  Kuchistu,  and  Y.  Morimo,  J.  Mol.  Spect. 

38,  70  (1971). 

i  R.  Walsh  and  S.  W.  Benson,  J.  Am.  Chem.  Soc.  88,  4570  (1966). 


FIGURE  1.  ILLUSTRATION  OF  THE  DEPENDENCE  OF  AN  AB  INITIO  CALCULATION 
ON  THE  BASIS  SET  AND  ON  THE  QUALITY  OF  THE  THEORY 
SZ,  SZP,  DZ,  AND  DZP  ARE  RESPECTIVELY  SINGLE  ZETA, 

SINGLE  ZETA  PLUS  POLARIZATION,  DOUBLE  ZETA,  ETC.  CONFIGURA¬ 
TION  INTERACTION  (Cl)  IS  USUALLY  ACCOMPLISHED  BY  ADDING 
SINGLE  AND  DOUBLE  EXCITATIONS.  MBPT  AND  CCM  IN  GENERAL 
EXCEED  SD-CI  IN  ACCURACY  SINCE  EFFECTS  OF  HIGHER  EXCITATIONS 
ARE  INCLUDED  TO  SOME  DEGREE.  MR-CCSD  INDICATES  COUPLED- 
CLUSTER  THEORY  LIMITED  TO  eT2+  T2  BUT  RELATIVE  TO  MORE  THAN 
ONE  REFERENCE  FUNCTION.  THE  BEST  POSSIBLE  SOLUTION  IN  A 
BASIS  SET  IS  FULL  Cl. 


AB  INITIO  CALCULATIONS 


\  Basis  Set  (Space) 
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ABSTRACT 


A  numerical  procedure  for  efficiently  solving  large  systems  of 
linear  equations  is  presented.  The  approach,  termed  the  reduced  linear 
equation  (RLE)  method,  is  illustrated  by  solving  the  systems  of  linear 
equations  that  arise  in  linearized  versions  of  coupled-cluster  theory. 
The  non-linear  coupled-cluster  equations  are  also  treated  with  the  RLE 
by  assuming  an  approximate  linearization  of  the  non-linear  terms.  Very 
efficient  convergence  for  linear  systems  and  good  convergence  for  non¬ 
linear  equations  is  found  for  a  numberof  examples  that  manifest  some 
degeneracy.  These  include  the  Be  atom,  H2  at  large  separation,  and  the 
N2  molecule.  The  RLE  method  is  compared  to  the  conventional  iterative 
procedure  and  to  Pade  approximants. 


I. 


INTRODUCTION 


The  basis  idea  of  the  coupled-cluster  method  (1-7)  relative  to 
a  single  reference  function  |$>,  is  that  the  exact  wavefunction  may  be 
written  as  f  =  exp(  T)|<j>>  where 

T-T,  ♦T2+...  (1) 


and  the  excitation  operators  {Tn }  are 


Tn  KKxc  -  Wi 


(2) 


By  back-projecting  y  onto  4>  and  a  sufficient  number  of  single,  double,  etc. 
excitations,  non-linear  algebraic  equations  of  the  form 


*  ?  Bijtj  *  ,/2!  CUkVk  +  1/31  ih  ■  0  (3) 

occur.  The  quantities  {t,-}  are  the  amplitudes  (tf *  *  * }  of  Eq.  (2)  which 

are  to  be  determined,  while  a.,  B. .,  C...,  etc.,  are  simply  combinations  of 

'  T  J  1 JK 

molecular  integrals.  Unlike  eigenvalue  equations,  the  coupled-cluster 
equations  are  independent  of  the  energy  attesting  to  the  linked-diagram, 
size-extensive  natureof  the  theory.^ 

In  most  of  the  applications  which  have  been  madif"10^he  trial  function 

^CCD  =  exP(T2^i>>  ^as  been  emPloyed,  which  terminates  Eq.  (3)  after  quadratic 
terms, 

'CCD  -  o+vvzi)!*"  .  (4) 

This  model  is  referred  to  as  coupled-cluster  doubles  (CCD)P’9^  or  by  Cizek  as 

(3  4) 

coupled-pair  many-electron  theory  (CPMET).  * 


tat 


2 


Another  model  called  coupled-cluster  singles  and  doubles  CCSD 

VT2 

can  be  defined  by  ^qqsd  =  e  I1*1*,  which  requires  that  T-j  be  included 
through  the  quartic  terms  while  T ^  still  appears  only  quadratically. 

Since  it  is  well-known  that  can  be  completely  eliminated  by  a  transforma¬ 
tion  to  Brueckner  orbitals  and  that  T]  is  usually  small  for  closed-shell 
(or  UHF  open-shell)  SCF  orbitals,  the  full  CCSD  model  appears  to  be 
unnecessarily  complex.  Consequently,  the  simpler  approximation  which 
includes  T-j  only  linearly,  and  which  we  will  refer  to  as  CCSD-1,  is 
considered.  (In  Paldus,  et.  al(^  ^  this  is  referred  to  as  approximation  B). 

The  trial  function  then  becomes, 

*CCS0-1  ■  (’  ♦  T1  ♦  T2  ♦  V2t!>I*»  •  (5) 

In  Eqs.  (4)  and  (5)  a  linearized  version  excluding  the  quadratic 

2 

terms  T 2  has  also  been  considered,  both  as  a  model  and  as  a  first  approximation 

to  the  solution  of  the  non-linear  equations.  In  the  case  of  CCD,  this  L-CCD 

model  corresponds  to  the  sum  of  all  double-excitation  diagrams  to  all  orders, 

while  L-CCSD  includes  the  single  excitation  diagram  as  well.  Both  models  are 

(11,12) 

size-extensive  as  required  by  the  linked-diagram  theorem.  If  the  L-CCD 

or  L-CCSD  wavefunction  were  used  in  an  expectation  value  formula  to  obtain  the 

amplitudes,  the  non-size-extensive  D-CI  and  SD-CI  models  would  then  emerge. 

2 

In  general,  the  quadratic  term  results  in  a  net  positive  contribution  as 
explained  elsewhere!  hence  the  L-CCD  and  L-CCSD  models,  though  size- 
extensive,  often  provide  correlation  energies  which  are  too  low  for  the 
basis  set.  Since  most  of  the  other  neglected  terms  like  Tj  also  result  in 
negative  values,  however,  the  net  error  is  normally  small,  except  for 
pathological  cases  usually  Involving  near-degenercies  where  the  damping 
effect  of  the  term  is  quite  important  7,13) 


3 


Perhaps  the  most  straightforward  approach  to  the  solution  of 

Eq.  (3)  is  to  employ  iterative  techniques.  Such  successive  iterations  of 

Eq.  (3)  lead  to  various  terms  in  the  linked-diagram  perturbation  expansion, 

( 3  1 

thereby  providing  a  means  to  sum  classes  of  MBPT  diagrams  to  all  orders.'  ' 
This  iterative  method  of  solution  makes  it  convenient  to  solve  the  coupled- 
cluster  equations  with  the  same  techniques  that  are  used  to  sum  MBPT 

(to) 

diagrams.  1,0  1  The  rate  of  convergence  can  usually  be  significantly 
enhanced  by  employing  Pade  approximantsl  M6ih  ich  are  simply  resummations 
of  the  energies  that  come  from  successive  iterations  of  the  equations. 

However,  there  can  be  some  difficulties  with  a  straight-forward  iterative 
approach. 

For  example,  if  the  reference  function  used  in  the  theory  were 

comparatively  poor  it  is  unreasonable  to  anticipate  that  good  convergence 

using  a  first-order  method  will  occur.  This  has  been  demonstrated  for  N2^7^ 

and  Be2^ 7  \  and  will  occur  in  almost  any  case  where  near  degeneracies  are 
(1 3 1 

a  problem.  '  A  common  example  in  molecular  theory  is  found  when  employing  a 

restricted  Hartree-Fock  (RHF)  solution  as  a  reference  function  at  large 

internuclear  separations,  where  the  RHF  function  erroneously  separates  to 

2 

ionic  products.  H2  is  a  good  example  of  this  problem,  since  the  lau 

2 

configuration  becomes  degenerate  with  the  la  RHF  configuration  at  large 

(18-20) 

separations.  Also,  Eq.  (3)  has  many  solutions',  so  convergence  to 
undesirable  solutions  is  also  possible.  In  some  cases  pertaining  to 
excited  states  and  excitation  energies,  it  is  necessary  to  obtain  these 
other  solutions  to  Eq.  (3)/19’^ 

An  alternative  method  of  solving  Eq.  (3)  for  the  CCD  case,  has  also 
i  5  ) 

been  used.  '  This  method  neglects  the  quadratic  term  initially,  solves 
the  set  of  linear  (L-CCD)  equations  exactly  (perhaps  by  matrix  inversion  or 


4 


by  iteration),  and  then  uses  a  Newton-Raphson  technique  for  the  quadratic 
(5) 

part.  '  A  few  repetitions  of  this  procedure  provides  the  solution.  The 
main  objection  to  this  approach  is  that  the  linear  terms  tend  to  be  negative 
while  the  non-linear  contributions  to  the  solution  are  positive.  A  simul¬ 
taneous  solution  embracing  both  terms  seems  to  be  preferable  to  exploit  the 
partial  cancellation.  In  particular,  an  initial  exact  solution  of  the  L-CCD 
equations  tends  to  provide  { t . }  that  are  often  far  away  from  the  correct  CCD 

J 

{tj}.  In  fact,  the  largest  differences  between  CCD  solutions  and  the  L-CCD 
solutions  occur  for  pathological  cases^’1'*’^  which  is  exactly  where  improved 
methods  of  solution  are  required. 

In  an  attempt  to  improve  the  convergence  of  the  solution  to  the 

coupled  cluster  equations,  we  have  investigated  a  technique  originally 

proposed  for  configuration  interaction  (Cl)  eigenvalue  equations  termed 

(21 -24) 

the  reduced  partitioning  (RP)  procedure.' 

In  the  Cl  case,  the  basic  idea  of  this  reduced  partitioning 

technique  is  to  drastically  reduce  the  effective  dimension  of  a  Cl  eigenvalue 

problem,  which  is  typically  very  large,  by  making  a  rectangular  transformation 

of  the  Cl  hamiltonian  matrix  to  a  set  of  m  trial  functions  (where  m  is  small). 

The  similarity  with  moment  theory[25^  the  Lanczos  algorithm!26^  and  Davidson's 

(27) 

method  for  eigenvalue  problems'  is  evident.  In  the  following  sections  we  show 
that  it  is  also  possible  to  use  a  reduce  partitioning  approach  for  linear 
equations.  We  also  show  that  this  reduced  linear  equation  method  can  be  used 
to  improve  convergence  for  the  non-linear  equations  of  the  type  that  occur  in 
coupled  cluster  theory. 

In  the  following,  this  technique  is  described  and  illustrated 
by  application  to  a  few  simple  cases. 
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IK _ THE  REDUCED  LINEAR  EQUATION  METHOD 

The  linear  approximation  to  the  coupled  cluster  equation  (3)  can 
be  written  in  matrix  form  as 


a  +  Bt  =  0 

O*  r\/\j  'l i 


(6) 


where  £  and  £  are  column  vectors  and  B  is  a  diagonally  dominant  square  matrix. 

In  the  L-CCD  model,  t  will  be  a  vector  with  a  length  equal  to  the  number  of 

double  excitations.  For  a  closed  shell  molecule  with  n  electrons  and  N  basis 
2  2 

functions,  roughly  n  N  excitations  are  possible.  Even  if  a  molecule  has 

twenty  electrons  and  one  hundred  basis  functions,  vectors  the  length  of  t  can 

be  conveniently  managed  on  most  medium-sized  computers. 

On  the  other  hand,  ^  is  a  square  matrix  and  the  number  of  elements 

which  must  be  processed  depends  upon  the  square  of  the  length  of  t.  In  the 

4  4 

L-CCD  model,  we  would  have  to  consider  up  to  n  N  elements  in  B  -  despite 

/V» 

the  fact  that  most  of  the  elements  are  zero.  For  the  model  problem  of  twenty 
electrons  and  a  hundred  basis  functions,  the  construction,  storage  and  manipula¬ 
tion  of  B  as  a  matrix  becomes  impractical. 

4 

Fortunately,  B  is  constructed  from  a  matrix  containing  only  N 
elements,  i.e.,  the  molecular  integrals,  and  efficient  programs  have  been 

4 

written  which  construct  the  products,  like  directly  from  the  N  elements. 

This  is  implicit  in  most  many-body  approaches  and  has  also  been  exploited 

(28  29)  4  4 

in  direct  □  techniques.  *  Thus,  the  n  N  bottleneck  is  broken  by  never 
constructing  B.  Instead,  a  program  taking  as  input  a  trial  vector  V  '  and 
the  integrals  directly  yields  a  new  quantity  which  is  used  to 

find  the  solution  Since  jj3  is  not  available,  most  standard  methods  for 
solving  systems  of  linear  equations  are  not  applicable. 
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One  simple  iterative  scheme  can  be  derived  by  partitioning  £  into 
a  diagonal  matrix  D  and  a  non-diagonal  matrix  Since  D  has  only  the  same 
number  of  elements  as  t,  D  can  also  be  managed  on  a  computer.  With  this  partition 

'Vi 

ing,  Ea.  (2)  can  be  rewritten  as 


et  *  -«  -«  t 
t  -  -e'1  k  +  iV 


Eq.  (4)  is  easily  changed  into  a  first-order  iterative  equation 

t<-i (»♦» tw) 

where  t^  is  the  m^  iteration  of  Eq.  (5)  and  where  usually  =  0. 

If  all  the  diagonal  terms  in  £  are  placed  into  J)  so  that  the 
diagonal  of  A  is  zero,  then  Eq.  (5)  is  the  Jacobi  method  for  solving 

linear  equations  IlnHoir  thoco  mnHitinnc  Fn.  1C  the  camo  ac 


Under  these  conditions,  Eq.  (9)  is  the  same  as 


the  perturbation  equation  for  the  L-CCD  model  assuming  an  Epstein-Nesbet 

fill 

partitioning  of  the  Hamiltonian.  J 1 ' 

(32) 

In  the  past,  the  M011er-Plesset  partitioning  has  been  demonstrated 
to  be  the  best  choice  for  solving  perturbation  or  cluster  equations.  The 
Mfiller-Plesset  partitioning  places  the  diagonal  effective  one-body  terms 
(e.g.,  orbital  energies)  into  D  and  leaves  the  two-body  terms  in  a.  The 
resulting  iterative  method  is  similar  to  the  Jacobi  method  with  over  relaxation, 
however,  here  there  is  a  different  relaxation  parameter  associated 
with  each  diagonal  element. 

The  slow  convergence  of  the  Jacobi  method  and  related  methods  is 


well-known.  To  obtain  more  rapid  covergence  for  the  iterative  solution  of  the 
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coupled  cluster  equations,  we  have  adapted  the  reduced  partitioning  (RP) 

(21-24) 

method  previously  applied  to  the  eigenvalue  problem.  the  RP  method, 

the  approximations  obtained  at  each  iteration  are  saved  and  then  used  on 
subsequent  iterations  to  provide  a  subpace  onto  which  Eq.  (6)  is 
projected  and  in  which  the  projected  system  of  equations  is  solved. 

For  example,  if  m  are  available,  then  we  set  up  a 

system  of  m-1  equations 


Rt  =  a 

Oj  'X/  'v 


(10) 


where 


r  =  y  t(i)  r  1 
ij  fa  k  Bk-e  t 


(j) 


-  i  ‘I1’  <°w  *  v  is) 


k  l 


(ID 


=  [  tk^  '  [  tk  “kk'k 


.(i) 


D,.f  t, 


(j+1) 


-  I  t 


(i) 


k  ak 


and  where 


“i  *  l  t 
1  k 


(i) 


k  ak 


.  ( i ) . . 


Once  t  is  determined,  the  best  approximation  to  x  using  m  V  '*s  is 


(12) 


z  • 

i-1 ,  m-1 


(13) 


The  superscript  [m]  on  t  indicates  that  x  has  been  obtained  from  m  t^'‘s. 

will  be  the  same  length  as  t.  As  can  be  seen  from  Eq.  (8)  and  (9),  the 
matrix  elements  for  R  and  a  are  easy  to  evaluate  since  they  are  only  weighted 


overlap  integrals  between  trial  solutions  of  different  iterations,  giving 
Eq.  (10.  This  small  matrix  problem 
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can  then  be  solved  with  standard  algorithms  specifically  designed  for 
approximately  singular  matrices.  It  is  important  to  note  that  the  R  matrix 
rapidly  tends  to  singularity.  However,  when  R  is  singular  then  t^  has 
probably  remained  unchanged  for  at  least  one  iteration  and  we  already  have 
a  solution. 

Notice  that  it  is  optional  whether  the  vectors  t^  are  not 

solutions  or  whether  increments  (i.e.,  successively  higher  orders  "perturbation" 

corrections  to  the  vector)  are  used,  since  the  two  approaches  are  related  by 

a  linear  transformation.  The  latter  approach  is  the  choice  in  the  RP 

(21,22) 

eigenvalue  problem  since  the  matrix  elements  in  Eq.  (10)  could  be 
conveniently  related  to  perturbation  energies  and  overlaps  of  perturbed  wave- 
functions. 

In  the  following  examples,  the  Jacobi  type  iteration  scheme  (c.f., 

Eq.  (9)  )has  been  used  to  generate  the  basis  {t^}  in  which  the  reduced 

'V  j 

linear  Eq.  (10)  is  solved.  The{£^}  could  provide  an  alternative  basis; 

j 

however,  the  convergence  properties  for  the  L-CCD  and  L-CCSD  models  would  remain 

j 

unaltered  and  a  time  consuminq  step-the  construction  of  {tLn-l}  would  have  been  added.  ! 

*\#  ■ 

Thus  far  we  have  described  methods  for  solving  linear  equations. 

The  general  coupled-cluster  equations  are  non-linear  and  it  is  of  interest 
to  ask  whether  the  reduced  linear  equation  method  can  be  used  for  solving  a 
matrix  equation  of  the  form  of  Eq.  (  3  ).  As  a  simple  first  step,  the 
quadratic  term  in  Eq.  (3)  is  written  as  an  effective  linear  term 


so  that  a  general  pseudo-linear  version  of  Eq.  (3)  limited  to  quadratic  terms 
can  be  v/rittcn  as 


a  +  B  (t)  t  =  0 

*  w  'X*  a. 


(15) 
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where  (5  (£)  is  |  +  ^(t).  If  we  partition  S  (t)  into  a  diagonal  matrix  J)  which 
does  not  depend  upon  t  and  a  non-diagonal  matrix  A ( t)  which  does  depend  upon 
t,  then  equations  (9)  ,(11 )  ,(12) ,  and  (13)  can  still  be  used  providing  A(t^) 
is  used  whenever  &  is  specified.  If  the  sequence 


( 16) 


converges,  then  the  sequence  of  solutions 

um.  i[2] . Im)  07) 

to  the  reduced  linear  equations 

R{m)(t)TCm]  *  £(m)  m=l,2,  3,  ...M  (18) 

will  converge  to  the  same  t  as  the  sequence  in  Eq.  (14).  If  the  iterative 
solution  to  Eq.  (9)  diverges,  then  the  sequence  of  solutions  to  the  reduced 
linear  equations  does  not  necessarily  converge  since  B  (t)  and  R(£)  do  not 
converge.  Consequently,  for  the  non-linear  equations,  it  is  useful  to  take 
the  best  available  solution  (i.e.,  t^)  and  use  it  to  define  B(t^).  Thus, 
a  suggested  process  for  non-linear  equations  is 

t(M+,)  =  <’[*  +  *  (*[m])  £(m)]  •  09) 

Alternative  schemes  for  handling  the  problems  of  non-linearity  in  a  linear 
framework  may  be  envisioned.  Such  schemes  must  address  questions  concerning 
consistency  of  the  final  solution  and  the  rate  of  convergence.  In  the  following, 
we  will  apply  the  simple  linearization  method  of  Eqs .(14)  and  (15)  to  several 
linear  and  non-linear  coupled  cluster  examples  to  illustrate  the  relative  speed 
of  convergence  using  the  reduced  linear  equation  (RLE)  approach  and  for  the 
simple  linearization  discussed  above  for  the  non-linear  equations. 
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III.  DISCUSSION  OF  RESULTS 

As  a  first  example  of  the  RLE  technique  consider  the  Be  atom. 

In  Table  I  are  shown  results  for  the  two  linearized  coupled-cluster  models, 

L-CCD  and  L-CCSD  and  the  non-linear  CCSD-1  model.  The  GOTO  7s3p  basis  is 

(33  7)  -5 

the  same  that  has  been  used  previously'  *  1  and  is  within  9x1  O  h  of  the 

2  2 

SCF  limit.  Since  in  Be,  the  2s  and  2p  orbitals  are  close  in  energy,  the  Is  2p 

2  2 

configuration  is  expected  to  be  relatively  important  compared  to  the  Is  2s 

reference  configuration  in  the  correlated  wavefunction.  In  the  intermediately 

2  2 

normalized  wavefunction,  the  correlating  Is  2p  configurations  have  coefficients 
of  0.14  for  CCSD-1  and  0.16  for  L-CCSD  while  other  double-excitation  configura¬ 
tions  involving  this  2p  orbital  have  coefficients  about  half  as  large.  The 
remaining  coefficients  are  <0.05.  This  might  be  compared  with  a  problem  like 

H20  at  equilibrium  where  the  largest  coefficient  of  any  configuration  is 
(81 

0.049.  This  feature  manifests  itself  in  comparatively  slow  convergence  of 
perturbation  theory  for  Be,  as  seen  in  Table  I.  On  the  other  hand,  the  Pade 
approximant  resummation  is  far  more  effective,  and  the  reduced  linear  equation 
approach  for  the  linearized  models  exceeds  the  convergence  of  the  Pade 
approximants. 

The  first  cycle  of  perturbation  theory  is  the  second-order  Mtfller- 
Plesset  energy^2)  while  at  least  the  third-order  perturbation  result  needs 
to  be  computed  before  either  the  [1,0]  Pade  approximant  [i.e.,  E2/(1-Ej/E2) , 
the  geometric  approximation]  or  the  first  cycle  of  the  RLE  is  possible.  In  fact, 
for  the  linear  theories  it  is  possible  to  get  two  perturbation  energies  per 
iteration,  using  a  version  of  the  2n+l  rule  of  perturbation  theoryP^  This  is 
not  in  general  applicable  to  the  non-linear  coupled-cluster  iterations,  however. 


-4 
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The  similarity  between  E2+E3>  obtained  from  the  second-cycle  of 
perturbation  theory,  and  the  [1,0]  approximant  attests  to  the  relative 
smallness  in  this  example  of  E3/E2  and  the  other  higher-order  differences 
that  distinguish  the  [1,0]  approximant  from  E2  +  E3<  Anologously,  the  Pade 
approximants  and  the  RLE  results  are  also  closely  related,  differing  primarily 
in  that  an  exact  matrix  problem  is  solved  in  the  RLE  case,  while  the  Pade 
approximants  offer  an  increasingly  good  approximation  to  the  matrix  equation 
solution. 

For  the  non-linear  CCSD-1  equations,  again  the  Pade  analysis  and 

the  RLE  analysis,  assuming  the  simple  linearization  of  Eq.  (14)  are  vastly 

superior  to  the  order-by-order  iterative  solution.  However,  the  first-order 

iterative  convergence  for  CCSD-1  is  superior  to  that  observed  in  the 

linearized  models.  This  is  because  the  incorporation  of  the  quadratic,  and 
2 

positive,  T2  term,  into  the  iterative  cycle  enhances  the  speed  of  convergence, 

2 

since  T2  acts  as  a  damping  factor  on  the  regative  linear  terms.  In  practice, 
two  iterations  of  the  linearized  equation  are  made,  followed  by  the  first  non¬ 
linear  iteration.  All  subsequent  cycles  embrace  a  linear  and  non-linear  iteration. 

The  main  advantage  of  the  Pade  approximant  analysis  compared  to 

the  RLE,  is  that  the  Pade  approximants  are  obtained  solely  from  the  energies 

(16) 

computed  at  each  iteration.  Once  the  energy  is  calculated,  previous  t  vectors 
can  be  discarded.  In  the  RLE,  the  matrix  in  Eq.  (10)  needs  to  be  computed, 
and  this  requires  overlaps  among  the  different  iterative  cycles,  which  requires 

a 

that  these  quantities  be  retained  on  mass  storage.  The  Pade  approximant 
analysis  also  has  its  disadvantages,  however,  in  that  they  do  not  provide 
wavefunction  information  conveniently,  i.i  the  RLE,  the  solution  of  the  matrix 
problem  gives  a  highly  accurate  representation  of  the  wavefunction  along  with 
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the  energy.  In  cases  where  t  has  not  yet  converged  the  RLE  solution  is 
available  and  exactly  corresponds  to  the  RLE  energy.  The  RLE  solution  is 
also  compact  in  the  sense  that  it  consists  of  just  a  few  coefficients  for 
the  different  t  vectors.  This  can  sometimes  have  advantages  when  using  the 
wavefunction  in  different  contexts,  such  as  the  prediction  of  second-order 
properties. ^ 

Another  difficult  example  for  convergence  of  the  coupled-cluster 
equations  is  offered  by  at  large  separation.  In  this  case,  the  reference 

p 

log  configuration  becomes  degenerate  with  the  correlating  double  excitation, 

p 

la£  .  Any  attempt  to  determine  such  a  potential  curve  with  only  a  single 
RHF  reference  function  and  finite  order  perturbation  theory  will  become 
suspect  at  sufficiently  large  internuclear  separation. 

This  is  illustrated  very  clearly  in  Tables  II  and  III,  where  results 
for  R  =  1.2  and  R  =  6.0  a.u.  are  compared.  The  L-CCD  results  converge  reasonably 
well  in  the  first  case  (-vlyh),  but  remain  0.4  hartrees  in  error  through  the 
same  cycle  at  6 a.u.  This  behavior  is  also  reflected  in  the  coefficient  of 
the  Ict^  configuration  which  is  0.06  at  1.2  a.u.  and  0.9  at  6.0  a.u.  It  is  also 
evident  that  the  L-CCD  result  is  much  too  low  at  6.0  a.u.  compared  to  the  CCD 
and  CCSD-1  values,  and  in  fact,  will  tend  to  minus  infinity  at  large  R  due 
to  singularities  in  the  energy  denominator.  On  the  other  hand,  the  L-CCD 
result  is  not  very  different  than  the  CCD  or  CCSD-1  values  at  1.2  a.u.  For  a 

problem  with  the  degeneracy  illustrated  here,  it  is  extremely  important  to 

2  (  7  1 

retain  the  quadratic  T2  terms  to  provide  a  reliable  answer.'  ' 

As  is  found  for  Be,  the  non-linear  CCD  and  CCSD-1  models  show  markedly 
better  convergence  than  t!,e  L-CCD  -model .  Whereas  at  1.2  a.u.  the  improvement  is 
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comparatively  slight,  at  6.0  a.u.,  there  is  a  dramatic  difference.  In  fact, 
the  first-order  iterative  CCD  solution  converges  as  fast  as  the  RLE  or  Pade 

approximants.  The  CCSD-1  has  somewhat  poorer  convergence  probably  due  to 

2 

the  failure  to  include  and  into  the  model. 

Figure  1  compares  the  size  of  the  increments  between  successive 
cycles  for  the  different  models,  the  first-order  iterative  approach,  the 
[N,N-1]  and  [N,N]  Pade  approximants,  and  the  RLE  results  for  the  L-CCD  and 
CCD  solutions  for  hL  at  R=6.0  a.u.  Rather  than  showing  the  convergence  to 
a  final  result,  this  figure  perhaps  offers  some  idea  of  the  stability  offered 
by  the  different  methods  used  to  converge  the  solutions.  The  solid  points 
refer  to  the  non-linear  CCD  results,  utilizing  the  simple  linearization 
techniques  of  Eq.  (14). 

It  is  clear  that  the  most  stable  results  even  for  this  highly 

degenerate  problem  are  offered  by  the  RLE  solution  to  the  L-CCD  equations. 

The  Pade  approximants  are  very  good,  but  they  still  tend  to  fluctuate  at  the 
-8 

+10  level.  The  standard  iterative,  perturbation  procedure  retains  large 

increments  at  about  the  10"1  level  in  all  cycles  shown. 

For  the  non-linear  CCD  equations,  there  is  little  difference 

between  the  RLE,  Pade,  and  straightforward  iterative  techniques.  None  of 

the  approaches  offer  as  much  numerical  stability  as  the  RLE  does  for  the 

linear  problem.  Rather  the  RLE  and  the  Pade  approximants  tend  to  reduce 
_g 

the  increment  to  ^10  in  the  8-10  cycles,  while  even  the  iterative  procedure 
-8 

goes  down  to  <10  .  Again  it  is  very  clear  that  the  iterative  approach  is 

vastly  superior  for  CCD  compared  to  L-CCD  for  this  pathological  example. 
Strictly  speaking,  the  [N,N-1]  and  [N,N]  Pade  approximants  correspond  to 
two  separate  sequences,  and  their  increments  should  be  compared  in  that 
•  manner,  but  for  this  example. there  is  no  appreciable  •difference . 
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Very  accurate  convergence  of  the  CCD  solution  has  been  found  to  be 
necessary  to  obtain  the  highly  sensitive  quartic  force  field  that  we  have 
reported  for  the  molecule/  '  In  this  example,  convergence  to  less 
than  ID'10  could  be  readily  detected  in  the  force  field  determination. 

In  a  final  example,  in  Table  IV,  results  for  the  CCD  and  CCSD-1 
models  for  N2  are  presented.  At  R=3.0  a.u.,  the  bifurcation  into  separate 
RHF  and  UHF  solutions  for  N2  has  already  occurred/17^  causing  a  decided 
lack  of  stability  in  the  RHF  based  coupled-cluster  models.  For  example,  some 
of  the  double  excitation  coefficients  are  already  as  large  as  0.22,  with 
several  others  lying  in  the  range  0.11  to  0.06. 

Except  for  the  much  larger  number  of  electrons  and  the  earlier 
onset  of  the  instability,  N2  behaves  very  much  like  H2,  in  that  an  L-CCD 
solution,  even  at  3.0  a.u.  will  be  much  too  negative  and  highly  slowly  convergent, 
while  the  CCD  and  CCSD-1  provide  more  realistic  results.  Again,  the  Pade 
approximants  and  RLE  give  relatively  good  convergence  to  the  solution,  both  better 
than  the  iterative  technique.  The  larger  number  of  electrons  does  not  seem  to 
require  that  more  cycles  be  run  to  get  the  solution,  even  though  it  is  clear 
that  the  iterative  technique  is  poorer  for  N2  than  in  H2. 

In  conclusion,  the  reduced  linear  equation  approach  has  been  shown 
to  be  very  efficient,  accurate,  and  convenient  for  obtaining  solutions  of 
linear  coupled  cluster  equations.  This  may  have  importance  for  solving  a 
wider  variety  of  linear  equation  systems.  The  RLE  coupled  with  a  simple, 
approximate  linearization  of  the  non-linear  equations,  is  somewhat  less 
satisfactory,  but  still  at  least  as  good  as  other  techniques,  and  much  better 
than  a  straightforward  iterative  approach.  More  sophisticated  ways  of 
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incorporating  the  non-linearity  of  the  general  coupled  cluster  equations 
that  still  give  a  quasi -linear  system  may  be  found  that  would  allow  the 
RLE  approach  to  function  more  efficiently.  The  main  advantage  of  the  RLE 
procedure  compared  to  a  Pade  approximant  analysis,  is  the  ease  of  obtaining 
the  coupled-cluster  wavefunction  in  addition  to  the  energy. 
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Figure  1.  Increments  between  successive  cycles  in  the  convergence  of 
the  L-CCD  and  CCD  equations  for  at  R  =  6.0  a.u.  Solid 
points  refer  to  the  L-CCD  equations  and  open  points  to  the 
CCD  equations. 


Open  points  are  CCD  and  solid 
points  are  L-CCD 
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